#### Code Preamble Read-Me ####

# This R Script contains the code needed to replicate all analysis in the forthcoming manuscript "Forecasting Partisan Collective Accountability During the 2024 U.S. Elections" to be published as a part of the Symposium on Forecasting the 2024 U.S. Elections in PS: Political Science & Politics.

# Note that this replication file assumes all source data files are saved in a folder called "Posted_Replication". Note all figures and tables will saved to this folder.

# The output of this R file is as follows and can be saved in the main replication directory:

# Manuscript:

#Figure 1: Presidential Approval & Incumbent Party Congressional Generic Percentage
#Figure 2: Marginal Effect of Presidential Approval & Party Brands on Election Outcomes
#Figure 3: Forecasting Model Out-of-Sample Predictions & Accuracy
#Table 1: 2024 Presidential Popular Vote Prediction Over Presidential Approval Levels
#Table 2: 2024 Presidential Electoral Vote Prediction Over Presidential Approval Level
#Table 3: 2024 U.S. Senate Prediction Over Generic Ballot Levels
#Table 4: 2024 U.S. House Prediction Over Generic Ballot Levels

# Appendix:

#Table A1: Quarterly Descriptive Statistics, 1937-2024
#Figure A1: Presidential Approval & Incumbent Party Congressional Generic Percentage
#Figure A2: Correlation of Quarterly Presidential Approval & Incumbent Party Brand By President
#Table A2: Quarterly OLS Models Predicting Incumbent Party Electoral Brand, 1937-2024
#Table A3: Annual Descriptive Statistics, 1938-2022
#Table A4: OLS Models Predicting Incumbent Party Popular Vote Percentage, 1940-2020
#Table A5: OLS Models Predicting Incumbent Party Electoral Votes, 1940-2020
#Table A6: OLS Models Predicting Incumbent Party U.S. Senate Seats, 1938-2022
#Table A7: OLS Models Predicting Incumbent Party U.S. House Seats, 1938-2022
#Figure A3: Marginal Effect of Presidential Approval on Outcomes By Election Type, 1940-2020
#Table A8: Out of Sample Predictions for Incumbent Party Popular Vote Percentage, 1940-2020
#Table A9: Out of Sample Predictions for Incumbent Party Electoral Votes, 1940-2020
#Table A10: Out of Sample Predictions for Incumbent Party U.S. Senate Seats, 1938-2022
#Table A11: Out of Sample Predictions for Incumbent Party U.S. House Seats, 1938-2022

#### Load Packages ####

library(readxl)
library(marginaleffects)
library(ggplot2)
library(modelsummary)
library(sandwich)
library(zoo)
library(DataCombine)
library(plyr)
library(broom)
library(lmtest)

#### Quarterly Data Management ####

load(file="/Users/caalgara/Desktop/dataverse_files/estimated_quarterly_generic_ballot.Rdata")

load(file="/Users/caalgara/Desktop/dataverse_files/estimated_presidential_approval.Rdata")

bootstrapped_dyad_ratios <- subset(bootstrapped_dyad_ratios,select=c(quarter,sd,lower,upper,bootstrapped_generic_estimate,date,year,decade)) #Note raw generic ballot differential is the quarterly Democratic-Republican differential. This must be recoded in the direction of the incumbent party.

colnames(bootstrapped_dyad_ratios) <- c("quarter","generic_ballot_sd","generic_ballot_lower","generic_ballot_upper","bootstrapped_generic_estimate","date","year","decade")

pres_bootstrapped_dyad_ratios <- subset(pres_bootstrapped_dyad_ratios,select=c(quarter,sd,lower,upper,bootstrapped_generic_estimate))
colnames(pres_bootstrapped_dyad_ratios) <- c("quarter","pres_approval_sd","pres_approval_lower","pres_approval_upper","bootstrapped_pres_approval")

gb_pres_dyad_ratios <- merge(bootstrapped_dyad_ratios,pres_bootstrapped_dyad_ratios,by=c("quarter"),all=T);rm(bootstrapped_dyad_ratios,pres_bootstrapped_dyad_ratios)

gb_pres_dyad_ratios$dem_incumbent_party <- ifelse(gb_pres_dyad_ratios$year %in% c(1938:1952,1961:1969,1977:1980,1993:2000,2009:2016,2021:2022),1,0)

gb_pres_dyad_ratios$incumbent_party_generic_ballot <- ifelse(gb_pres_dyad_ratios$dem_incumbent_party == 1,gb_pres_dyad_ratios$bootstrapped_generic_estimate,ifelse(gb_pres_dyad_ratios$dem_incumbent_party==0,100-gb_pres_dyad_ratios$bootstrapped_generic_estimate,NA))

cor(gb_pres_dyad_ratios$incumbent_party_generic_ballot,gb_pres_dyad_ratios$bootstrapped_pres_approval,use="complete.obs")

gb_pres_dyad_ratios$pres_administration <- ifelse(gb_pres_dyad_ratios$year %in% c(1937:1944),"Franklin Delano Roosevelt, 1937-1945",ifelse(gb_pres_dyad_ratios$quarter %in% c(1945.1),"Franklin Delano Roosevelt, 1937-1945",ifelse(gb_pres_dyad_ratios$quarter %in% c(1945.2,1945.3,1945.4),"Harry S. Truman, 1945-1952",ifelse(gb_pres_dyad_ratios$year %in% c(1946:1952),"Harry S. Truman, 1945-1952",ifelse(gb_pres_dyad_ratios$year %in% c(1953:1960),"Dwight Eisenhower, 1953-1960",ifelse(gb_pres_dyad_ratios$year %in% c(1961:1963),"John F. Kennedy, 1961-1963",ifelse(gb_pres_dyad_ratios$year %in% c(1964:1968),"Lyndon B. Johnson, 1963-1968",ifelse(gb_pres_dyad_ratios$year %in% c(1969:1973),"Richard Nixon, 1969-1974",ifelse(gb_pres_dyad_ratios$quarter %in% c(1974.1,1974.2,1974.3),"Richard Nixon, 1969-1974",ifelse(gb_pres_dyad_ratios$quarter %in% c(1974.4),"Gerald Ford, 1974-1976",ifelse(gb_pres_dyad_ratios$year %in% c(1975:1976),"Gerald Ford, 1974-1976",ifelse(gb_pres_dyad_ratios$year %in% c(1977:1980),"Jimmy Carter, 1977-1980",ifelse(gb_pres_dyad_ratios$year %in% c(1981:1988),"Ronald Reagan 1981-1988",ifelse(gb_pres_dyad_ratios$year %in% c(1989:1992),"George H.W. Bush, 1989-1992",ifelse(gb_pres_dyad_ratios$year %in% c(1993:2000),"Bill Clinton, 1993-2000",ifelse(gb_pres_dyad_ratios$year %in% c(2001:2008),"George W. Bush, 2001-2008",ifelse(gb_pres_dyad_ratios$year %in% c(2009:2016),"Barack Obama, 2009-2016",ifelse(gb_pres_dyad_ratios$year %in% c(2017:2020),"Donald Trump, 2017-2020","Joe Biden, 2021-2024"))))))))))))))))))

gb_pres_dyad_ratios$pres_administration <- factor(gb_pres_dyad_ratios$pres_administration,levels=c("Franklin Delano Roosevelt, 1937-1945","Harry S. Truman, 1945-1952","Dwight Eisenhower, 1953-1960","John F. Kennedy, 1961-1963","Lyndon B. Johnson, 1963-1968","Richard Nixon, 1969-1974","Gerald Ford, 1974-1976","Jimmy Carter, 1977-1980","Ronald Reagan 1981-1988","George H.W. Bush, 1989-1992","Bill Clinton, 1993-2000","George W. Bush, 2001-2008","Barack Obama, 2009-2016","Donald Trump, 2017-2020","Joe Biden, 2021-2024"))

r_df <- ddply(gb_pres_dyad_ratios, .(pres_administration), summarise, rsq=round(cor(bootstrapped_pres_approval,incumbent_party_generic_ballot,use="complete.obs"),2))

plot <- ggplot(gb_pres_dyad_ratios,aes(x=bootstrapped_pres_approval,y=incumbent_party_generic_ballot,label=quarter)) + geom_point(shape=21) + theme_bw() + facet_wrap(~pres_administration,scales="free") + stat_smooth(method="lm") + scale_x_continuous("Quarterly Presidential Approval",limits=c(34,70),breaks=seq(35,70,5),labels=paste(seq(35,70,5),"%",sep="")) + scale_y_continuous("Quarterly Incumbent Party Generic Ballot Electoral Brand",limits=c(34,70),breaks=seq(35,70,5),labels=paste(seq(35,70,5),"%",sep=""))+ geom_text(data=r_df, aes(label=paste("corr=", rsq)),x=-Inf, y=Inf, hjust=-0.2, vjust=1.2);print(plot)
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/FigA2.png", width = 12, height = 8, units = "in",bg="white")

r_df <- ddply(gb_pres_dyad_ratios, .(), summarise, rsq=round(cor(bootstrapped_pres_approval,incumbent_party_generic_ballot,use="complete.obs"),3))

model1 <- lm(incumbent_party_generic_ballot~bootstrapped_pres_approval,data=gb_pres_dyad_ratios)
x <- coeftest(model1,vcov=vcovHC(model1,"HC2"))       
tidy(x,conf.int=T) # coefficient, R^2, & standard error reported in caption of Figure 1 
prediction <- predictions(model1,vcov="HC2")

plot <- ggplot(gb_pres_dyad_ratios,aes(x=bootstrapped_pres_approval,y=incumbent_party_generic_ballot)) + geom_point(shape=21) + theme_bw() + scale_x_continuous("Quarterly Presidential Approval",limits=c(34,70),breaks=seq(35,70,5),labels=paste(seq(35,70,5),"%",sep="")) + scale_y_continuous("Quarterly Incumbent Party Generic \nBallot Electoral Brand",limits=c(34,70),breaks=seq(35,70,5),labels=paste(seq(35,70,5),"%",sep="")) + geom_text(data=r_df, aes(label=paste("corr=", rsq)),x=-Inf, y=Inf, hjust=-0.2, vjust=1.2) + geom_line(data=prediction,aes(x=bootstrapped_pres_approval,y=estimate),color="blue",inherit.aes = F) + geom_ribbon(data=prediction,aes(x=bootstrapped_pres_approval,y=estimate,ymin=conf.low,ymax=conf.high),alpha=0.2)
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/Fig1.png", width = 6, height = 3.25, units = "in",bg="white") 

gb_pres_dyad_ratios$year_fe <- factor(gb_pres_dyad_ratios$year)

plot <- ggplot(gb_pres_dyad_ratios,aes(x=quarter,y=bootstrapped_pres_approval)) + geom_line() + theme_bw() + scale_x_continuous("",breaks=seq(1940,2024,10)) + scale_y_continuous("U.S. Presidential Approval",breaks=seq(-100,100,5),labels=paste(seq(-100,100,5),"%",sep=""),limits=c(34,70)) + geom_hline(yintercept=50,color="red",linetype="dashed")
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/FigA1a.png", width = 8, height = 4, units = "in",bg="white")

gb_pres_dyad_ratios$incumbent_party_generic_ballot_lower <- ifelse(gb_pres_dyad_ratios$dem_incumbent_party == 1,gb_pres_dyad_ratios$generic_ballot_lower,ifelse(gb_pres_dyad_ratios$dem_incumbent_party==0,100-gb_pres_dyad_ratios$generic_ballot_lower,NA))

gb_pres_dyad_ratios$incumbent_party_generic_ballot_upper <- ifelse(gb_pres_dyad_ratios$dem_incumbent_party == 1,gb_pres_dyad_ratios$generic_ballot_upper,ifelse(gb_pres_dyad_ratios$dem_incumbent_party==0,100-gb_pres_dyad_ratios$generic_ballot_upper,NA))

plot <- ggplot(gb_pres_dyad_ratios,aes(x=quarter,y=incumbent_party_generic_ballot)) + geom_line() + theme_bw() + scale_x_continuous("",breaks=seq(1940,2024,10)) + scale_y_continuous("Incumbent Party Generic Ballot Percentage",breaks=seq(-100,100,2),labels=paste(seq(-100,100,2),"%",sep=""),limits=c(40,60)) + geom_hline(yintercept=50,color="red",linetype="dashed")
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/FigA1b.png", width = 8, height = 4, units = "in",bg="white")

#### Merge with Quarterly Covariate Data ####

gdp <- read.csv("/Users/caalgara/Desktop/dataverse_files/GDPC1-2.csv")

gdp_annual <- read.csv("/Users/caalgara/Desktop/dataverse_files/gdp_growth_annual.csv")

unemployment <- read.csv("/Users/caalgara/Desktop/dataverse_files/unemployment.csv")

time_change <- read.csv("/Users/caalgara/Desktop/dataverse_files/Incumbent party in power time.csv")
time_change <- na.omit(time_change)

colnames(unemployment)[1] <- c("quarter")

econ <- merge(unemployment,time_change,by="quarter",all=T)

gdp$quarter <- as.character(as.yearqtr(gdp$DATE, format = "%Y-%m-%d"))
gdp$quarter  <- gsub("Q","",gdp$quarter)
gdp$quarter  <- gsub(" ",".",gdp$quarter)
gdp$year <-  substr(gdp$quarter,1,4)
gdp$quarterly_gdp <- gdp$GDPC1

gdp <- subset(gdp,select=c(quarter,year,quarterly_gdp))
gdp_annual$DATE <- substr(gdp_annual$DATE,1,4)
colnames(gdp_annual) <- c("year","annual_gdp_growth")

gdp <- slide(gdp, Var = "quarterly_gdp", slideBy = -1,NewVar = "lagged_quarterly_gdp")
gdp$quarterly_gdp_growth <- (gdp$quarterly_gdp-gdp$lagged_quarterly_gdp)/gdp$quarterly_gdp

econ <- merge(econ,gdp,by=c("quarter"),all=T)
econ$year <- substr(econ$quarter,1,4)

econ <- merge(econ,gdp_annual,by="year",all=T)
econ <- subset(econ,!is.na(econ$quarter))

rm(time_change,unemployment,gdp,gdp_annual)

econ$year <- NULL

gb_pres_dyad_ratios <- merge(gb_pres_dyad_ratios,econ,by=c("quarter"),all=T)

######### Party Brands Presidential Approval Models #########

summary(model1 <- lm(incumbent_party_generic_ballot ~ bootstrapped_pres_approval + time_party_in_power_presidency + year_fe,data=gb_pres_dyad_ratios))

summary(model2 <- lm(incumbent_party_generic_ballot ~ bootstrapped_pres_approval + time_party_in_power_presidency + pres_administration,data=gb_pres_dyad_ratios))

summary(model3 <- lm(incumbent_party_generic_ballot ~ bootstrapped_pres_approval + time_party_in_power_presidency + quarterly_gdp_growth + unemployment_rate + year_fe,data=gb_pres_dyad_ratios))

summary(model4 <- lm(incumbent_party_generic_ballot ~ bootstrapped_pres_approval + time_party_in_power_presidency + quarterly_gdp_growth + unemployment_rate + pres_administration,data=gb_pres_dyad_ratios))

models <- list(model1,model2,model3,model4)

cm <- c('bootstrapped_pres_approval'='Presidential Approval','time_party_in_power_presidency'='Time in Power','quarterly_gdp_growth'='GDP Growth','unemployment_rate'='Unemployment Rate')

ar <- data.frame("FEs & Clustered SEs","Yearly","Administration","Yearly","Administration")

modelsummary(models,vcov = list(vcovCL(model1, cluster = ~year_fe),vcovCL(model2,cluster=~pres_administration),vcovCL(model3, cluster = ~year_fe),vcovCL(model4,cluster=~pres_administration)),gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",coef_map = cm,add_rows = ar,title="Quarterly OLS Models Predicting Incumbent Party Electoral Brand, 1937-2022",output = "/Users/caalgara/Desktop/dataverse_files/TableA2.html",notes="DV: Incumbent party electoral generic ballot brand.")

modelsummary(models,vcov = list(vcovCL(model1, cluster = ~year_fe),vcovCL(model2,cluster=~pres_administration),vcovCL(model3, cluster = ~year_fe),vcovCL(model4,cluster=~pres_administration)),gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",coef_map = cm,add_rows = ar,title="Quarterly OLS Models Predicting Incumbent Party Electoral Brand, 1937-2022",output = "/Users/caalgara/Desktop/dataverse_files/TableA2.tex",notes="DV: Incumbent party electoral generic ballot brand.")

x <- subset(gb_pres_dyad_ratios,select=c(incumbent_party_generic_ballot,bootstrapped_pres_approval,time_party_in_power_presidency,quarterly_gdp_growth,unemployment_rate))
colnames(x) <- c("Incumbent Party Party Brand","Presidential Approval","Quarterly Incumbent Party Time in Power","Quarterly GDP Growth","Quarterly Unemployment Rate")

datasummary(`Incumbent Party Party Brand` + `Presidential Approval` + `Quarterly Incumbent Party Time in Power` + `Quarterly GDP Growth` + `Quarterly Unemployment Rate` ~ N + Mean + Min + Max + Median + SD,data=x,title="Quarterly Descriptive Statistics, 1937-2022",fmt=2,na.rm=T,output="/Users/caalgara/Desktop/dataverse_files/TableA1.tex",notes = "GDP missing = 1937.3-1947.1 & 2024.3 (40 Qs) | Unemployment missing = 1937.3-1937.4 & 2024.3 (3 Qs)")

datasummary(`Incumbent Party Party Brand` + `Presidential Approval` + `Quarterly Incumbent Party Time in Power` + `Quarterly GDP Growth` + `Quarterly Unemployment Rate` ~ N + Mean + Min + Max + Median + SD,data=x,title="Quarterly Descriptive Statistics, 1937-2022",fmt=2,na.rm=T,output="/Users/caalgara/Desktop/dataverse_files/TableA1.html",notes = "GDP missing = 1937.3-1947.1 & 2024.3 (40 Qs) | Unemployment missing = 1937.3-1937.4 & 2024.3 (3 Qs)")


rm(ar,model1,model2,model3,model4,models,cm,plot,x,r_df,prediction)

#### Construct Election Datasets ####

pres <- read_excel("/Users/caalgara/Desktop/dataverse_files/forecasting_outcome_var_data_collection_pres.xlsx")
pres <- subset(pres,pres$year != 1936)

senate <- read_excel("/Users/caalgara/Desktop/dataverse_files/forecasting_outcome_var_data_collection_senate.xlsx",sheet = "senate")
colnames(senate) <- c("year","senate_dem_seats_won","senate_gop_seats_won","senate_dem_seats_before","senate_gop_seats_before")

house <- read.csv("/Users/caalgara/Desktop/dataverse_files/forecasting_outcome_var_data_collection_house.csv")
colnames(house) <- c("year","house_dem_seats_won","house_gop_seats_won","house_dem_seats_before","house_gop_seats_before")

returns <- merge(senate,house,by=c("year"))
returns <- merge(pres,returns,by=c("year"),all=T)

# Taking measures third quarter of election year

x <- subset(gb_pres_dyad_ratios,gb_pres_dyad_ratios$quarter %in% seq(1938.3,2022.3,2))

returns <- merge(returns,x,by=c("year"))

returns$dem_pres_twoparty_vote <- (returns$dem_votes/(returns$dem_votes+returns$gop_votes)) * 100
returns$dem_incumbent_party <- factor(returns$dem_incumbent_party,labels=c("Republican","Democratic"))

returns$open_presidential_election <- ifelse(returns$year %in% c(1952,1960,1968,1988,2000,2008,2016),1,ifelse(returns$year %in% c(1940,1944,1948,1956,1964,1972,1976,1980,1984,1992,1996,2004,2012,2020),0,NA))

returns$incumbent_party_generic_ballot <- ifelse(returns$dem_incumbent_party == "Democratic",returns$bootstrapped_generic_estimate,100-returns$bootstrapped_generic_estimate)

returns$incumbent_party_twoparty_vote <- ifelse(returns$dem_incumbent_party == "Democratic",returns$dem_pres_twoparty_vote,100-returns$dem_pres_twoparty_vote)

returns$incumbent_party_ev <- ifelse(returns$dem_incumbent_party == "Democratic",returns$dem_electoral_votes,returns$gop_electoral_votes)

returns$incumbent_party_senate <- ifelse(returns$dem_incumbent_party == "Democratic",returns$senate_dem_seats_won,returns$senate_gop_seats_won)

returns$incumbent_party_house <- ifelse(returns$dem_incumbent_party=="Democratic",returns$house_dem_seats_won,returns$house_gop_seats_won)

returns$midterm_election <- ifelse(returns$year %in% seq(1938,2022,4),1,0)

returns$incumbent_party_house_seats_before <- ifelse(returns$dem_incumbent_party == "Democratic",returns$house_dem_seats_before,ifelse(returns$dem_incumbent_party == "Republican",returns$house_gop_seats_before,NA))

returns$incumbent_party_senate_seats_before <- ifelse(returns$dem_incumbent_party == "Democratic",returns$senate_dem_seats_before,ifelse(returns$dem_incumbent_party == "Republican",returns$senate_gop_seats_before,NA))

returns$incumbent_party_senate_seat_change <- returns$incumbent_party_senate-returns$incumbent_party_senate_seats_before

returns$incumbent_party_house_seat_change <- returns$incumbent_party_house-returns$incumbent_party_house_seats_before

returns$dem_incumbent_party <- factor(returns$dem_incumbent_party,levels=c("Republican","Democratic"))
returns$dem_incumbent_party_num <- as.numeric(returns$dem_incumbent_party)

rm(x)

#### Incumbent Presidential Election Models: OLS ####

summary(model2 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot,data=returns))
summary(model1 <- lm(incumbent_party_twoparty_vote~bootstrapped_pres_approval,data=returns))
summary(model3 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot+bootstrapped_pres_approval,data=returns))
summary(model4 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot+bootstrapped_pres_approval + dem_incumbent_party_num + time_party_in_power_presidency,data=returns))
summary(model5 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot+bootstrapped_pres_approval + dem_incumbent_party_num + time_party_in_power_presidency + unemployment_rate + annual_gdp_growth,data=returns))

models <- list(model1,model2,model3,model4,model5)

cm <- c('bootstrapped_pres_approval'='Presidential Approval','incumbent_party_generic_ballot'="Incumbent Party Brand",'dem_incumbent_party_num'="Democratic Administration",'time_party_in_power_presidency'='Time in Power','annual_gdp_growth'='GDP Growth','unemployment_rate'='Unemployment Rate')

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party Popular Vote Percentage, 1940-2020",coef_map = cm,notes = "DV: Two-party presidential vote-share won by incumbent party.",output = "/Users/caalgara/Desktop/dataverse_files/TableA4.html")

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party Popular Vote Percentage, 1940-2020",coef_map = cm,notes = "DV: Two-party presidential vote-share won by incumbent party.",output = "/Users/caalgara/Desktop/dataverse_files/TableA4.tex")

mes5_pop_vote <- marginaleffects::avg_slopes(model5,variables = c("incumbent_party_generic_ballot","bootstrapped_pres_approval"),vcov = "HC2")

# Electoral Votes

summary(model2 <- lm(incumbent_party_ev~incumbent_party_generic_ballot,data=returns))
summary(model1 <- lm(incumbent_party_ev~bootstrapped_pres_approval,data=returns))
summary(model3 <- lm(incumbent_party_ev~incumbent_party_generic_ballot+bootstrapped_pres_approval,data=returns))
summary(model4 <- lm(incumbent_party_ev~incumbent_party_generic_ballot+bootstrapped_pres_approval + dem_incumbent_party_num + time_party_in_power_presidency,data=returns))
summary(model5 <- lm(incumbent_party_ev~incumbent_party_generic_ballot+bootstrapped_pres_approval + dem_incumbent_party_num + time_party_in_power_presidency + unemployment_rate + annual_gdp_growth,data=returns))

models <- list(model1,model2,model3,model4,model5)

cm <- c('bootstrapped_pres_approval'='Presidential Approval','incumbent_party_generic_ballot'="Incumbent Party Brand",'dem_incumbent_party_num'="Democratic Administration",'time_party_in_power_presidency'='Time in Power','annual_gdp_growth'='GDP Growth','unemployment_rate'='Unemployment Rate')

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party Electoral Votes, 1940-2020",coef_map = cm,notes = "DV: Electoral votes won by incumbent party.","/Users/caalgara/Desktop/dataverse_files/TableA5.html")

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party Electoral Votes, 1940-2020",coef_map = cm,notes = "DV: Electoral votes won by incumbent party.","/Users/caalgara/Desktop/dataverse_files/TableA5.tex")

mes5_evs <- marginaleffects::avg_slopes(model5,variables = c("incumbent_party_generic_ballot","bootstrapped_pres_approval"),vcov = "HC2")

#### Heterogenity in Presidential Approval By Election Type ####

returns$open_presidential_election <- ifelse(returns$year %in% c(1952,1960,1968,1988,2000,2008,2016),1,ifelse(returns$year %in% c(1940,1944,1948,1956,1964,1972,1976,1980,1984,1992,1996,2004,2012,2020),0,NA))

summary(model1 <- lm(incumbent_party_twoparty_vote~bootstrapped_pres_approval*open_presidential_election,data=returns))

me1 <- marginaleffects::avg_slopes(model1,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

summary(model3 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot+bootstrapped_pres_approval*open_presidential_election,data=returns))

me3 <- marginaleffects::avg_slopes(model3,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

summary(model4 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot+bootstrapped_pres_approval*open_presidential_election + dem_incumbent_party_num + time_party_in_power_presidency,data=returns))

me4 <- marginaleffects::avg_slopes(model4,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

summary(model5 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot+bootstrapped_pres_approval*open_presidential_election + dem_incumbent_party_num + time_party_in_power_presidency + unemployment_rate + annual_gdp_growth,data=returns))

me5 <- marginaleffects::avg_slopes(model5,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

me1$model <- 1
me3$model <- 3
me4$model <- 4
me5$model <- 5

x <- rbind(me1,me3,me4,me5)

x$open_presidential_election <- factor(x$open_presidential_election,levels=c(0,1),labels=c("Incumbent \n Re-Election","Open \nSeat"))

x$model <- factor(x$model,levels=c(1,3,4,5),labels=c(expression(atop("Model 1 Specification", paste("Adjusted ",R^2," = 0.671"))),expression(atop("Model 3 Specification", paste("Adjusted ",R^2," = 0.657"))),expression(atop("Model 4 Specification", paste("Adjusted ",R^2," = 0.697"))),expression(atop("Model 5 Specification", paste("Adjusted ",R^2," = 0.664")))))

x$label <- ifelse(x$p.value < 0.10 & x$p.value > 0.05,paste(round(x$estimate,2),"*",sep=""),ifelse(x$p.value < 0.05 & x$p.value > 0.01,paste(round(x$estimate,2),"**",sep=""),ifelse(x$p.value < 0.01,paste(round(x$estimate,2),"***",sep=""),NA)))
x$label[x$label %in% "0.6***"] <- "0.60***"
x$label[x$label %in% "0.3***"] <- "0.30***"

plot <- ggplot(x,aes(x=open_presidential_election,y=estimate,label=label)) + geom_pointrange(aes(x= open_presidential_election, ymin = conf.low, ymax = conf.high), lwd = 0.75,fill="white",shape=21) + theme_bw() + coord_flip() + geom_text(vjust=-1,hjust=0.25) + geom_hline(yintercept = 0, colour = "red", lty = 2) + scale_x_discrete("") + scale_y_continuous("Average Marginal Effect of Presidential Approval on Incumbent Party Popular Vote-Percentage") + labs(caption="* p < 0.10 ** p < 0.05; *** p < 0.01. N = 21 Presidential Elections") + theme(axis.text.x = element_text(hjust = 0.5),axis.text.y = element_text(hjust = 0.5)) + facet_wrap(~model,scales="free",labeller = label_parsed);print(plot)
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/FigA3a.png", width = 8, height = 4, units = "in",bg="white")

# Electoral Votes

summary(model1 <- lm(incumbent_party_ev~bootstrapped_pres_approval*open_presidential_election,data=returns))

me1 <- marginaleffects::avg_slopes(model1,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

summary(model3 <- lm(incumbent_party_ev~incumbent_party_generic_ballot+bootstrapped_pres_approval*open_presidential_election,data=returns))

me3 <- marginaleffects::avg_slopes(model3,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

summary(model4 <- lm(incumbent_party_ev~incumbent_party_generic_ballot+bootstrapped_pres_approval*open_presidential_election + dem_incumbent_party_num + time_party_in_power_presidency,data=returns))

me4 <- marginaleffects::avg_slopes(model4,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

summary(model5 <- lm(incumbent_party_ev~incumbent_party_generic_ballot+bootstrapped_pres_approval*open_presidential_election + dem_incumbent_party_num + time_party_in_power_presidency + unemployment_rate + annual_gdp_growth,data=returns))

me5 <- marginaleffects::avg_slopes(model5,variables = c("bootstrapped_pres_approval"),by="open_presidential_election",vcov = "HC2")

me1$model <- 1
me3$model <- 3
me4$model <- 4
me5$model <- 5

x <- rbind(me1,me3,me4,me5)

x$open_presidential_election <- factor(x$open_presidential_election,levels=c(0,1),labels=c("Incumbent \n Re-Election","Open \nSeat"))

x$model <- factor(x$model,levels=c(1,3,4,5),labels=c(expression(atop("Model 1 Specification", paste("Adjusted ",R^2," = 0.606"))),expression(atop("Model 3 Specification", paste("Adjusted ",R^2," = 0.586"))),expression(atop("Model 4 Specification", paste("Adjusted ",R^2," = 0.582"))),expression(atop("Model 5 Specification", paste("Adjusted ",R^2," = 0.539")))))

x$label <- ifelse(x$p.value < 0.10 & x$p.value > 0.05,paste(round(x$estimate,2),"*",sep=""),ifelse(x$p.value < 0.05 & x$p.value > 0.01,paste(round(x$estimate,2),"**",sep=""),ifelse(x$p.value < 0.01,paste(round(x$estimate,2),"***",sep=""),NA)))
x$label[x$label %in% "7.7*"] <- "7.70*"

plot <- ggplot(x,aes(x=open_presidential_election,y=estimate,label=label)) + geom_pointrange(aes(x= open_presidential_election, ymin = conf.low, ymax = conf.high), lwd = 0.75,fill="white",shape=21) + theme_bw() + coord_flip() + geom_text(vjust=-1,hjust=0.25) + geom_hline(yintercept = 0, colour = "red", lty = 2) + scale_x_discrete("") + scale_y_continuous("Average Marginal Effect of Presidential Approval on Incumbent Party Popular Vote-Percentage") + labs(caption="* p < 0.10 ** p < 0.05; *** p < 0.01. N = 21 Presidential Elections") + theme(axis.text.x = element_text(hjust = 0.5),axis.text.y = element_text(hjust = 0.5)) + facet_wrap(~model,scales="free",labeller = label_parsed);print(plot)
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/FigA3b.png", width = 8, height = 4, units = "in",bg="white")

#### Incumbent U.S. Senate Election Models: OLS ####

summary(model2 <- lm(incumbent_party_senate~incumbent_party_generic_ballot,data=returns))
summary(model1 <- lm(incumbent_party_senate~bootstrapped_pres_approval,data=returns))
summary(model3 <- lm(incumbent_party_senate~incumbent_party_generic_ballot+bootstrapped_pres_approval,data=returns))
summary(model4 <- lm(incumbent_party_senate~incumbent_party_generic_ballot+bootstrapped_pres_approval + midterm_election,data=returns))
summary(model5 <- lm(incumbent_party_senate~incumbent_party_generic_ballot+bootstrapped_pres_approval + midterm_election + dem_incumbent_party_num + time_party_in_power_presidency,data=returns))
summary(model6 <- lm(incumbent_party_senate~incumbent_party_generic_ballot+bootstrapped_pres_approval + midterm_election + dem_incumbent_party_num + time_party_in_power_presidency + unemployment_rate + annual_gdp_growth,data=returns))

models <- list(model1,model2,model3,model4,model5,model6)

cm <- c('bootstrapped_pres_approval'='Presidential Approval','incumbent_party_generic_ballot'="Incumbent Party Brand",'midterm_election'='Midterm Election','dem_incumbent_party_num'="Democratic Administration",'time_party_in_power_presidency'='Time in Power','annual_gdp_growth'='GDP Growth','unemployment_rate'='Unemployment Rate')

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party U.S. Senate Seats, 1938-2022",coef_map = cm,notes = "DV: U.S. Senate seats won by incumbent party.","/Users/caalgara/Desktop/dataverse_files/TableA6.html")

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party U.S. Senate Seats, 1938-2022",coef_map = cm,notes = "DV: U.S. Senate seats won by incumbent party.","/Users/caalgara/Desktop/dataverse_files/TableA6.tex")

mes6_senate <- marginaleffects::avg_slopes(model6,variables = c("incumbent_party_generic_ballot","bootstrapped_pres_approval"),vcov = "HC2")

#### Incumbent U.S. House Election Models: OLS ####

summary(model2 <- lm(incumbent_party_house~incumbent_party_generic_ballot,data=returns))
summary(model1 <- lm(incumbent_party_house~bootstrapped_pres_approval,data=returns))
summary(model3 <- lm(incumbent_party_house~incumbent_party_generic_ballot+bootstrapped_pres_approval,data=returns))
summary(model4 <- lm(incumbent_party_house~incumbent_party_generic_ballot+bootstrapped_pres_approval + midterm_election,data=returns))
summary(model5 <- lm(incumbent_party_house~incumbent_party_generic_ballot+bootstrapped_pres_approval + midterm_election + dem_incumbent_party_num + time_party_in_power_presidency,data=returns))
summary(model6 <- lm(incumbent_party_house~incumbent_party_generic_ballot+bootstrapped_pres_approval + midterm_election + dem_incumbent_party_num + time_party_in_power_presidency + unemployment_rate + annual_gdp_growth,data=returns))

models <- list(model1,model2,model3,model4,model5,model6)

cm <- c('bootstrapped_pres_approval'='Presidential Approval','incumbent_party_generic_ballot'="Incumbent Party Brand",'midterm_election'='Midterm Election','dem_incumbent_party_num'="Democratic Administration",'time_party_in_power_presidency'='Time in Power','annual_gdp_growth'='GDP Growth','unemployment_rate'='Unemployment Rate')

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party U.S. House Seats, 1938-2022",coef_map = cm,notes = "DV: U.S. House seats won by incumbent party.","/Users/caalgara/Desktop/dataverse_files/TableA7.html")

modelsummary(models,vcov = "HC2",gof_omit = 'DF|Deviance|RMSE|AIC|BIC|F|Log.Lik|Std.Errors',fmt = 2,estimate = "{estimate}{stars}",title="OLS Models Predicting Incumbent Party U.S. House Seats, 1938-2022",coef_map = cm,notes = "DV: U.S. House seats won by incumbent party.","/Users/caalgara/Desktop/dataverse_files/TableA7.tex")

mes6_house <- marginaleffects::avg_slopes(model6,variables = c("incumbent_party_generic_ballot","bootstrapped_pres_approval"),vcov = "HC2")

#### OLS Plots of Models & Descriptive Statistics ####

mes5_pop_vote$model <- "Presidential Popular Vote"
mes5_evs$model <- "Presidential Electoral Votes"
mes6_senate$model <- "U.S. Senate Seats"
mes6_house$model <- "U.S. House Seats"

x <- rbind(mes5_pop_vote,mes5_evs,mes6_senate,mes6_house)

x$label <- ifelse(x$p.value < 0.10 & x$p.value > 0.05,paste(round(x$estimate,2),"*",sep=""),ifelse(x$p.value < 0.05 & x$p.value > 0.01,paste(round(x$estimate,2),"**",sep=""),ifelse(x$p.value < 0.01,paste(round(x$estimate,2),"***",sep=""),NA)))
x$label[x$label == "1.5***"] <- "1.50***"
x$model <- factor(x$model,levels=c("Presidential Popular Vote","Presidential Electoral Votes","U.S. Senate Seats","U.S. House Seats"))

x$var <- ifelse(x$term %in% "bootstrapped_pres_approval","Presidential Job \nApproval",ifelse(x$term %in% "incumbent_party_generic_ballot","Incumbent Party \nBrand",NA))
x$var <- factor(x$var,levels=c("Incumbent Party \nBrand","Presidential Job \nApproval"))

plot <- ggplot(x,aes(x=var,y=estimate,label=label)) + geom_pointrange(aes(x= var, ymin = conf.low, ymax = conf.high), lwd = 0.75,fill="white",shape=21) + theme_bw() + coord_flip() + geom_text(vjust=-1,hjust=0.25) + geom_hline(yintercept = 0, colour = "red", lty = 2) + scale_x_discrete("") + scale_y_continuous("Average Marginal Effect of Incumbent Party Evaluations on Outcomes") + labs(caption="* p < 0.10 ** p < 0.05; *** p < 0.01.") + theme(axis.text.x = element_text(hjust = 0.5),axis.text.y = element_text(hjust = 0.5)) + facet_wrap(~model,scales="free");print(plot)
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/Fig2.png", width = 8, height = 4, units = "in",bg="white")

x <- subset(returns,select=c(incumbent_party_twoparty_vote,incumbent_party_ev,incumbent_party_senate,incumbent_party_house,incumbent_party_generic_ballot,bootstrapped_pres_approval,midterm_election,dem_incumbent_party_num,time_party_in_power_presidency,annual_gdp_growth,unemployment_rate))

colnames(x) <- c("Incumbent Party Two-Party Vote % Won","Incumbent Party Electoral Votes Won","Incumbent Party U.S. Senate Seats Won","Incumbent Party U.S. House Seats Won","Incumbent Party Party Brand","Presidential Approval","Midterm Election Cycle","Democratic Administration","Quarterly Incumbent Party Time in Power","Quarterly GDP Growth","Annual Unemployment Rate")

datasummary(`Incumbent Party Two-Party Vote % Won` + `Incumbent Party Electoral Votes Won` + `Incumbent Party U.S. Senate Seats Won` + `Incumbent Party U.S. House Seats Won` + `Incumbent Party Party Brand` + `Presidential Approval` + `Midterm Election Cycle` + `Democratic Administration` + `Quarterly Incumbent Party Time in Power` + `Quarterly GDP Growth` + `Annual Unemployment Rate` ~ N + Mean + Min + Max + Median + SD,data=x,title="Annual Descriptive Statistics, 1938-2022",fmt=2,na.rm=T,output="/Users/caalgara/Desktop/dataverse_files/TableA3.html")

datasummary(`Incumbent Party Two-Party Vote % Won` + `Incumbent Party Electoral Votes Won` + `Incumbent Party U.S. Senate Seats Won` + `Incumbent Party U.S. House Seats Won` + `Incumbent Party Party Brand` + `Presidential Approval` + `Midterm Election Cycle` + `Democratic Administration` + `Quarterly Incumbent Party Time in Power` + `Quarterly GDP Growth` + `Annual Unemployment Rate` ~ N + Mean + Min + Max + Median + SD,data=x,title="Annual Descriptive Statistics, 1938-2022",fmt=2,na.rm=T,output="/Users/caalgara/Desktop/dataverse_files/TableA3.tex")

#### Incumbent OLS Election Models: Out of Sample Predictions ####

out_of_sample_predictions <- list()

for(i in seq(1940,2020,4)){
  y <- subset(returns,returns$year == i)
  x <- subset(returns,returns$year != i)
  
  summary(model4 <- lm(incumbent_party_twoparty_vote~incumbent_party_generic_ballot+bootstrapped_pres_approval,data=x))
  
  os <- data.frame(year=i,adjusted_r2 = summary(model4)$adj.r.squared)
  
  prediction <- avg_predictions(model4,vcov="HC2", newdata = datagrid(incumbent_party_generic_ballot = y$incumbent_party_generic_ballot, bootstrapped_pres_approval = y$bootstrapped_pres_approval, grid_type = "counterfactual"))
  
  os$prediction <- prediction$estimate
  os$prediction_lower <- prediction$conf.low
  os$prediction_upper <- prediction$conf.high
  
  os$observed_result <- y$incumbent_party_twoparty_vote
  os$error <- os$observed_result - os$prediction
  os$absolute_error <- abs(os$error)
  
  os$correctly_predicted <- ifelse(os$observed_result >= 50 & os$prediction >= 50, "Yes",ifelse(os$observed_result < 50 & os$prediction < 50, "Yes","No"))
  
  out_of_sample_predictions[[i]] <- os
}

out_of_sample_predictions_pres_pop_vote <- ldply(out_of_sample_predictions,data.frame)

out_of_sample_predictions <- list()

for(i in seq(1940,2020,4)){
  y <- subset(returns,returns$year == i)
  x <- subset(returns,returns$year != i)
  
  summary(model1 <- lm(incumbent_party_ev~incumbent_party_generic_ballot+bootstrapped_pres_approval,data=x))
  
  os <- data.frame(year=i,adjusted_r2 = summary(model1)$adj.r.squared)
  
  prediction <- avg_predictions(model1,vcov="HC2", newdata = datagrid(incumbent_party_generic_ballot = y$incumbent_party_generic_ballot, bootstrapped_pres_approval = y$bootstrapped_pres_approval, grid_type = "counterfactual"))
  
  os$prediction <- prediction$estimate
  os$prediction_lower <- prediction$conf.low
  os$prediction_upper <- prediction$conf.high
  
  os$observed_result <- y$incumbent_party_ev
  os$error <- os$observed_result - os$prediction
  os$absolute_error <- abs(os$error)
  
  os$correctly_predicted <- ifelse(os$observed_result >= 270 & os$prediction >= 270, "Yes",ifelse(os$observed_result < 270 & os$prediction < 270, "Yes","No"))
  
  out_of_sample_predictions[[i]] <- os
}

out_of_sample_predictions_pres_ev <- ldply(out_of_sample_predictions,data.frame)

out_of_sample_predictions <- list()

for(i in seq(1938,2022,2)){
  y <- subset(returns,returns$year == i)
  x <- subset(returns,returns$year != i)
  
  summary(model1 <- lm(incumbent_party_senate~incumbent_party_generic_ballot+bootstrapped_pres_approval+midterm_election,data=x))
  
  os <- data.frame(year=i,adjusted_r2 = summary(model1)$adj.r.squared)
  
  prediction <- avg_predictions(model1,vcov="HC2", newdata = datagrid(incumbent_party_generic_ballot = y$incumbent_party_generic_ballot, bootstrapped_pres_approval = y$bootstrapped_pres_approval,midterm_election=y$midterm_election, grid_type = "counterfactual"))
  
  os$prediction <- prediction$estimate
  os$prediction_lower <- prediction$conf.low
  os$prediction_upper <- prediction$conf.high
  
  os$observed_result <- y$incumbent_party_senate
  os$error <- os$observed_result - os$prediction
  os$absolute_error <- abs(os$error)
  
  os$correctly_predicted <- ifelse(os$year < 1960 & os$observed_result >= 48 & os$prediction >= 48,"Yes",ifelse(os$year < 1960 & os$observed_result < 48 & os$prediction < 48,"Yes",ifelse(os$year >= 1960 & os$observed_result >= 50 & os$prediction >= 50,"Yes",ifelse(os$year >= 1960 & os$observed_result < 50 & os$prediction < 50,"Yes","No"))))
  
  out_of_sample_predictions[[i]] <- os
}

out_of_sample_predictions_senate <- ldply(out_of_sample_predictions,data.frame)

out_of_sample_predictions <- list()

for(i in seq(1938,2022,2)){
  y <- subset(returns,returns$year == i)
  x <- subset(returns,returns$year != i)
  
  summary(model1 <- lm(incumbent_party_house~incumbent_party_generic_ballot+bootstrapped_pres_approval+midterm_election,data=x))
  
  os <- data.frame(year=i,adjusted_r2 = summary(model1)$adj.r.squared)
  
  prediction <- avg_predictions(model1,vcov="HC2", newdata = datagrid(incumbent_party_generic_ballot = y$incumbent_party_generic_ballot, bootstrapped_pres_approval = y$bootstrapped_pres_approval,midterm_election=y$midterm_election, grid_type = "counterfactual"))
  
  os$prediction <- prediction$estimate
  os$prediction_lower <- prediction$conf.low
  os$prediction_upper <- prediction$conf.high
  
  os$observed_result <- y$incumbent_party_house
  os$error <- os$observed_result - os$prediction
  os$absolute_error <- abs(os$error)
  
  os$correctly_predicted <- ifelse(os$observed_result >= 218 & os$prediction >= 218, "Yes",ifelse(os$observed_result < 218 & os$prediction < 218, "Yes","No"))
  
  out_of_sample_predictions[[i]] <- os
}

out_of_sample_predictions_house <- ldply(out_of_sample_predictions,data.frame)

# Make Scattplots of Out-of-Sample-Predictions

out_of_sample_predictions_pres_pop_vote$model <- "Presidential Popular Vote \nAbsolute Error: Mean = 2.53, Median = 1.68, \nMin = 0.02, Max = 7.35, SD = 2.28 \nDiscrete Elections Predicted Correctly: 19/21"
out_of_sample_predictions_pres_ev$model <- "Presidential Electoral Votes \nAbsolute Error: Mean = 83.98, Median = 75.16, \nMin = 5.37, Max = 226.01, SD = 64.43 \nDiscrete Elections Predicted Correctly: 15/21"
out_of_sample_predictions_senate$model <- "U.S. Senate Seats \nAbsolute Error: Mean = 4.48, Median = 5.45, \nMin = 0.08, Max = 18.16, SD = 3.80 \nDiscrete Elections Predicted Correctly: 25/43"
out_of_sample_predictions_house$model <- "U.S. House Seats \nAbsolute Error: Mean = 17.40, Median = 17.11, \nMin =0.02, Max = 54.71,  SD = 12.83 \nDiscrete Elections Predicted Correctly: 35/43"

plot <- ggplot(out_of_sample_predictions_pres_pop_vote,aes(x=prediction,y=observed_result,label=year)) + geom_text() + theme_bw() + facet_wrap(~model,scales="free") + scale_x_continuous("Incumbent-Party Model Out-of-Sample Prediction",limits=c(42,62),breaks=seq(42,62,2),labels=paste(seq(42,62,2),"%",sep="")) + scale_y_continuous("Incumbent-Party Observed Result",limits=c(42,62),breaks=seq(42,62,2),labels=paste(seq(42,62,2),"%",sep="")) + geom_abline(intercept = 0, slope = 1) + annotate("text",x=45,y=60,label="italic(Over-Performance)",color="blue",parse=T) + annotate("text",x=60,y=45,label="italic(Under-\nPerformance)",color="red",parse=T) + labs(caption="Predictions derived from Model (3).")
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/Fig3A.png", width = 7, height = 4, units = "in",bg="white")

plot <- ggplot(out_of_sample_predictions_pres_ev,aes(x=prediction,y=observed_result,label=year)) + geom_text() + theme_bw() + facet_wrap(~model,scales="free") + scale_x_continuous("Incumbent-Party Model Out-of-Sample Prediction",limits=c(45,550),breaks=seq(50,550,50)) + scale_y_continuous("Incumbent-Party Observed Result",limits=c(45,550),breaks=seq(50,550,50)) + geom_abline(intercept = 0, slope = 1) + annotate("text",x=100,y=500,label="italic(Over-Performance)",color="blue",parse=T) + annotate("text",x=500,y=100,label="italic(Under-\nPerformance)",color="red",parse=T) + labs(caption="Predictions derived from Model (3).")
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/Fig3B.png", width = 7, height = 4, units = "in",bg="white")

plot <- ggplot(out_of_sample_predictions_senate,aes(x=prediction,y=observed_result,label=year)) + geom_text() + theme_bw() + facet_wrap(~model,scales="free") + scale_x_continuous("Incumbent-Party Model Out-of-Sample Prediction",limits=c(30,75),breaks=seq(30,75,5)) + scale_y_continuous("Incumbent-Party Observed Result",limits=c(30,75),breaks=seq(30,75,5)) + geom_abline(intercept = 0, slope = 1) + annotate("text",x=35,y=70,label="italic(Over-Performance)",color="blue",parse=T) + annotate("text",x=70,y=35,label="italic(Under-\nPerformance)",color="red",parse=T) + labs(caption="Predictions derived from Model (3).")
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/Fig3C.png", width = 7, height = 4, units = "in",bg="white")

plot <- ggplot(out_of_sample_predictions_house,aes(x=prediction,y=observed_result,label=year)) + geom_text() + theme_bw() + facet_wrap(~model,scales="free") + scale_x_continuous("Incumbent-Party Model Out-of-Sample Prediction",limits=c(130,300),breaks=seq(150,300,25)) + scale_y_continuous("Incumbent-Party Observed Result",limits=c(130,300),breaks=seq(150,300,25)) + geom_abline(intercept = 0, slope = 1) + annotate("text",x=150,y=280,label="italic(Over-Performance)",color="blue",parse=T) + annotate("text",x=280,y=150,label="italic(Under-\nPerformance)",color="red",parse=T) + labs(caption="Predictions derived from Model (3).")
ggsave(plot,file="/Users/caalgara/Desktop/dataverse_files/Fig3D.png", width = 7, height = 4, units = "in",bg="white")

# Tables of Out-of-Sample Predictions

out_of_sample_predictions_pres_pop_vote$model <- NULL
out_of_sample_predictions_pres_pop_vote$year <- as.character(out_of_sample_predictions_pres_pop_vote$year)

colnames(out_of_sample_predictions_pres_pop_vote) <- c("Election","Adjusted R2","Prediction","95% Lower Bound CI","95% Upper Bound CI","Observed Result","Residual Error","Absolute Error","Discrete Correct Prediction")
out_of_sample_predictions_pres_pop_vote <- out_of_sample_predictions_pres_pop_vote[,c(1,2,6,3,4,5,7,8,9)]
datasummary_df(out_of_sample_predictions_pres_pop_vote,title="Out of Sample Predictions for Incumbent Party Popular Vote Percentage, 1940-2020",output = "/Users/caalgara/Desktop/dataverse_files/TableA8.html",notes="DV: Two-party presidential vote-share won by incumbent party. 19/21 correctly predicted. Predictions derived from Model (3).")

datasummary_df(out_of_sample_predictions_pres_pop_vote,title="Out of Sample Predictions for Incumbent Party Popular Vote Percentage, 1940-2020",output = "/Users/caalgara/Desktop/dataverse_files/TableA8.tex",notes="DV: Two-party presidential vote-share won by incumbent party. 19/21 correctly predicted. Predictions derived from Model (3).")

out_of_sample_predictions_pres_ev$model <- NULL
out_of_sample_predictions_pres_ev$year <- as.character(out_of_sample_predictions_pres_ev$year)

colnames(out_of_sample_predictions_pres_ev) <- c("Election","Adjusted R2","Prediction","95% Lower Bound CI","95% Upper Bound CI","Observed Result","Residual Error","Absolute Error","Discrete Correct Prediction")
out_of_sample_predictions_pres_ev <- out_of_sample_predictions_pres_ev[,c(1,2,6,3,4,5,7,8,9)]
datasummary_df(out_of_sample_predictions_pres_ev,title="Out of Sample Predictions for Incumbent Party Electoral Votes, 1940-2020",output = "/Users/caalgara/Desktop/dataverse_files/TableA9.html",notes="DV: Electoral votes won by incumbent party. 15/21 correctly predicted. Predictions derived from Model (3).")

datasummary_df(out_of_sample_predictions_pres_ev,title="Out of Sample Predictions for Incumbent Party Electoral Votes, 1940-2020",output = "/Users/caalgara/Desktop/dataverse_files/TableA9.tex",notes="DV: Electoral votes won by incumbent party. 15/21 correctly predicted. Predictions derived from Model (3).")

out_of_sample_predictions_senate$model <- NULL
out_of_sample_predictions_senate$year <- as.character(out_of_sample_predictions_senate$year)

colnames(out_of_sample_predictions_senate) <- c("Election","Adjusted R2","Prediction","95% Lower Bound CI","95% Upper Bound CI","Observed Result","Residual Error","Absolute Error","Discrete Correct Prediction")
out_of_sample_predictions_senate <- out_of_sample_predictions_senate[,c(1,2,6,3,4,5,7,8,9)]
datasummary_df(out_of_sample_predictions_senate,title="Out of Sample Predictions for Incumbent Party U.S. Senate Seats, 1938-2022",output = "/Users/caalgara/Desktop/dataverse_files/TableA10.html",notes="DV: U.S. Senate seats won by incumbent party. 25/43 correctly predicted. Predictions derived from Model (4).")

datasummary_df(out_of_sample_predictions_senate,title="Out of Sample Predictions for Incumbent Party U.S. Senate Seats, 1938-2022",output = "/Users/caalgara/Desktop/dataverse_files/TableA10.tex",notes="DV: U.S. Senate seats won by incumbent party. 25/43 correctly predicted. Predictions derived from Model (4).")

out_of_sample_predictions_house$model <- NULL
out_of_sample_predictions_house$year <- as.character(out_of_sample_predictions_house$year)

colnames(out_of_sample_predictions_house) <- c("Election","Adjusted R2","Prediction","95% Lower Bound CI","95% Upper Bound CI","Observed Result","Residual Error","Absolute Error","Discrete Correct Prediction")
out_of_sample_predictions_house <- out_of_sample_predictions_house[,c(1,2,6,3,4,5,7,8,9)]
datasummary_df(out_of_sample_predictions_house,title="Out of Sample Predictions for Incumbent Party U.S. House Seats, 1938-2022",output = "/Users/caalgara/Desktop/dataverse_files/TableA11.html",notes="DV: U.S. House seats won by incumbent party. 35/43 correctly predicted. Predictions derived from Model 43).")

datasummary_df(out_of_sample_predictions_house,title="Out of Sample Predictions for Incumbent Party U.S. House Seats, 1938-2022",output = "/Users/caalgara/Desktop/dataverse_files/TableA11.tex",notes="DV: U.S. House seats won by incumbent party. 35/43 correctly predicted. Predictions derived from Model (4).")

#### Incumbent Election Models: 2024 Predictions as of August 19, 2024 ####

predictions_2022_pop <- list()

for(i in seq(38,55,1)){
  
  summary(model_ols <- lm(incumbent_party_twoparty_vote~bootstrapped_pres_approval+ incumbent_party_generic_ballot,data=returns))
  
  prediction_ols <- avg_predictions(model_ols, vcov="HC2",newdata = datagrid(bootstrapped_pres_approval = i,incumbent_party_generic_ballot= (46.3/(46.3+45.2))*100,grid_type = "counterfactual"))
  prediction_ols$i <- i
  
  prediction <- data.frame(presidential_approval = i, ols_estimate = prediction_ols$estimate,ols_lb = prediction_ols$conf.low , ols_ub = prediction_ols$conf.high)
  
  predictions_2022_pop[[i]] <- prediction
}

predictions_2022_pop <- ldply(predictions_2022_pop,data.frame)

predictions_2022_ev <- list()

for(i in seq(38,55,1)){
  y <- returns
  
  summary(model_ols <- lm(incumbent_party_ev~bootstrapped_pres_approval+incumbent_party_generic_ballot,data=returns))
  
  prediction_ols <- avg_predictions(model_ols, vcov="HC2",newdata = datagrid(bootstrapped_pres_approval = i,incumbent_party_generic_ballot= (46.3/(46.3+45.2))*100,grid_type = "counterfactual"))
  
  prediction <- data.frame(presidential_approval = i, ols_estimate = prediction_ols$estimate,ols_lb = prediction_ols$conf.low , ols_ub = prediction_ols$conf.high)
  
  predictions_2022_ev[[i]] <- prediction
}

predictions_2022_ev <- ldply(predictions_2022_ev,data.frame)

#
predictions_2022_senate <- list()

for(i in seq(47,53,1)){
  
  summary(ols_model <- lm(incumbent_party_senate~incumbent_party_generic_ballot+bootstrapped_pres_approval+midterm_election,data=returns))
  
  prediction_ols <- avg_predictions(ols_model,vcov="HC2",newdata = datagrid(bootstrapped_pres_approval = (38.2/(38.2+55.8))*100,incumbent_party_generic_ballot= i,midterm_election=0,grid_type = "counterfactual"))
  
  prediction <- data.frame(generic_ballot = i, ols_estimate = prediction_ols$estimate,ols_lb = prediction_ols$conf.low , ols_ub = prediction_ols$conf.high)
  
  predictions_2022_senate[[i]] <- prediction
}

predictions_2022_senate <- ldply(predictions_2022_senate,data.frame)

#
predictions_2022_house <- list()

for(i in seq(47,53,1)){
  

  summary(ols_model <- lm(incumbent_party_house~incumbent_party_generic_ballot+bootstrapped_pres_approval+midterm_election,data=returns))
  
  prediction_ols <- avg_predictions(ols_model,vcov="HC2", newdata = datagrid(bootstrapped_pres_approval = (38.2/(38.2+55.8))*100,incumbent_party_generic_ballot= i,midterm_election=0,grid_type = "counterfactual"))
  
  prediction <- data.frame(generic_ballot = i, ols_estimate = prediction_ols$estimate,ols_lb = prediction_ols$conf.low , ols_ub = prediction_ols$conf.high)
  
  predictions_2022_house[[i]] <- prediction
}

predictions_2022_house <- ldply(predictions_2022_house,data.frame)

colnames(predictions_2022_pop) <- c("Presidential Approval Level","Percentage Estimate","95% Percentage Lower Bound CI","95% Percentage Upper Bound CI")

predictions_2022_pop <- predictions_2022_pop[,c(1:4)]
datasummary_df(predictions_2022_pop,title="2024 Presidential Popular Vote Prediction Over Presidential Approval Levels",output = "/Users/caalgara/Desktop/dataverse_files/Table1.html",notes="Predictions derived from Model (3) & observed covariate values on 8/19/2024.")

datasummary_df(predictions_2022_pop,title="2024 Presidential Popular Vote Prediction Over Presidential Approval Levels",output = "/Users/caalgara/Desktop/dataverse_files/Table1.tex",notes="Predictions derived from Model (3) & observed covariate values on 8/19/2024.")

colnames(predictions_2022_ev) <- c("Presidential Approval Level","Electoral Votes Won Estimate","95% Votes Lower Bound CI","95% Votes Upper Bound CI")

predictions_2022_ev <- predictions_2022_ev[,c(1:4)]
datasummary_df(predictions_2022_ev,title="2024 Presidential Electoral Vote Prediction Over Presidential Approval Level",output = "/Users/caalgara/Desktop/dataverse_files/Table2.html",notes="Predictions derived from Model (3) & observed covariate values on 8/19/2024.")

datasummary_df(predictions_2022_ev,title="2024 Presidential Electoral Vote Prediction Over Presidential Approval Level",output = "/Users/caalgara/Desktop/dataverse_files/Table2.tex",notes="Predictions derived from Model (3) & observed covariate values on 8/19/2024.")

colnames(predictions_2022_senate) <- c("Generic Ballot Support Level","Seats Won Estimate","95% Seats Lower Bound CI","95% Seats Upper Bound CI")

predictions_2022_senate <- predictions_2022_senate[,c(1:4)]
datasummary_df(predictions_2022_senate,title="2024 U.S. Senate Prediction Over Generic Ballot Levels",output = "/Users/caalgara/Desktop/dataverse_files/Table3.html",notes="Predictions derived from Model (4) & observed covariate values on 8/19/2024.")

datasummary_df(predictions_2022_senate,title="2024 U.S. Senate Prediction Over Generic Ballot Levels",output = "/Users/caalgara/Desktop/dataverse_files/Table3.tex",notes="Predictions derived from Model (4) & observed covariate values on 8/15/2024.")

colnames(predictions_2022_house) <-  c("Generic Ballot Support Level","Seats Won Estimate","95% Seats Lower Bound CI","95% Seats Upper Bound CI")

predictions_2022_house <- predictions_2022_house[,c(1:4)]
datasummary_df(predictions_2022_house,title="2024 U.S. House Prediction Over Generic Ballot Levels",output = "/Users/caalgara/Desktop/dataverse_files/Table4.html",notes="Predictions derived from Model (4) & observed covariate values on 8/15/2024.")

datasummary_df(predictions_2022_house,title="2024 U.S. House Prediction Over Generic Ballot Levels",output = "/Users/caalgara/Desktop/dataverse_files/Table4.tex",notes="Predictions derived from Model (4) & observed covariate values on 8/15/2024.")

